Давайте проанализируем данные опроса 4361 женщин из Ботсваны:
botswana.tsv О каждой из них мы знаем:
сколько детей она родила (признак ceb) возраст (age) длительность получения образования (educ) религиозная принадлежность (religion) идеальное, по её мнению, количество детей в семье (idlnchld) была ли она когда-нибудь замужем (evermarr) возраст первого замужества (agefm) длительность получения образования мужем (heduc) знает ли она о методах контрацепции (knowmeth) использует ли она методы контрацепции (usemeth) живёт ли она в городе (urban) есть ли у неё электричество, радио, телевизор и велосипед (electric, radio, tv, bicycle) Давайте научимся оценивать количество детей ceb по остальным признакам.
Загрузите данные и внимательно изучите их. Сколько разных значений принимает признак religion?
Во многих признаках есть пропущенные значения. Сколько объектов из 4361 останется, если выбросить все, содержащие пропуски?
В разных признаках пропуски возникают по разным причинам и должны обрабатываться по-разному.
Например, в признаке agefm пропуски стоят только там, где evermarr=0, то есть, они соответствуют женщинам, никогда не выходившим замуж. Таким образом, для этого признака NaN соответствует значению "не применимо".
В подобных случаях, когда признак x1 на части объектов в принципе не может принимать никакие значения, рекомендуется поступать так:
создать новый бинарный признак x2={1,0,x1='не применимо',иначе; заменить "не применимо" в x1 на произвольную константу c, которая среди других значений x1 не встречается. Теперь, когда мы построим регрессию на оба признака и получим модель вида y=β0+β1x1+β2x2, на тех объектах, где x1 было измерено, регрессионное уравнение примет вид y=β0+β1x, а там, где x1 было "не применимо", получится y=β0+β1c+β2. Выбор c влияет только на значение и интерпретацию β2, но не β1.
Давайте используем этот метод для обработки пропусков в agefm и heduc.
Создайте признак nevermarr, равный единице там, где в agefm пропуски. Удалите признак evermarr — в сумме с nevermarr он даёт константу, значит, в нашей матрице X будет мультиколлинеарность. Замените NaN в признаке agefm на cagefm=0. У объектов, где nevermarr = 1, замените NaN в признаке heduc на cheduc1=−1 (ноль использовать нельзя, так как он уже встречается у некоторых объектов выборки). Сколько осталось пропущенных значений в признаке heduc?
Избавимся от оставшихся пропусков.
Для признаков idlnchld, heduc и usemeth проведите операцию, аналогичную предыдущей: создайте индикаторы пропусков по этим признакам (idlnchld_noans, heduc_noans, usemeth_noans), замените пропуски на нехарактерные значения (cidlnchld=−1, cheduc2=−2 (значение -1 мы уже использовали), cusemeth=−1).
Остались только пропуски в признаках knowmeth, electric, radio, tv и bicycle. Их очень мало, так что удалите объекты, на которых их значения пропущены.
Какого размера теперь наша матрица данных? Умножьте количество строк на количество всех столбцов (включая отклик ceb).
Постройте регрессию количества детей ceb на все имеющиеся признаки методом smf.ols, как в разобранном до этого примере. Какой получился коэффициент детерминации R2? Округлите до трёх знаков после десятичной точки.
Проверьте критерием Бройша-Пагана гомоскедастичность ошибки в построенной модели. Выполняется ли она? Если ошибка гетероскедастична, перенастройте модель, сделав поправку Уайта типа HC1.
Удалите из модели незначимые признаки religion, radio и tv. Проверьте гомоскедастичность ошибки, при необходимости сделайте поправку Уайта.
Не произошло ли значимого ухудшения модели после удаления этой группы признаков? Проверьте с помощью критерия Фишера. Чему равен его достигаемый уровень значимости? Округлите до четырёх цифр после десятичной точки.
Если достигаемый уровень значимости получился маленький, верните все удалённые признаки; если он достаточно велик, оставьте модель без религии, тв и радио.
Признак usemeth_noans значим по критерию Стьюдента, то есть, при его удалении модель значимо ухудшится. Но вообще-то отдельно его удалять нельзя: из-за того, что мы перекодировали пропуски в usemeth произвольно выбранным значением cusemeth=−1, удалять usemeth_noans и usemeth можно только вместе.
Удалите из текущей модели usemeth_noans и usemeth. Проверьте критерием Фишера гипотезу о том, что качество модели не ухудшилось. Введите номер первой значащей цифры в достигаемом уровне значимости (например, если вы получили 5.5×10−8, нужно ввести 8).
Если достигаемый уровень значимости получился маленький, верните удалённые признаки; если он достаточно велик, оставьте модель без usemeth и usemeth_noans.
In [1]:
import statsmodels
import scipy as sc
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.stats.api as sms
from statsmodels.graphics.regressionplots import plot_leverage_resid2
import matplotlib.pyplot as plt
In [2]:
%pylab inline
In [49]:
raw = pd.read_csv("botswana.tsv", sep="\t", index_col=False)
raw.head()
Out[49]:
In [50]:
raw.religion.value_counts()
Out[50]:
In [51]:
raw.shape
Out[51]:
In [52]:
raw.dropna().shape
Out[52]:
In [53]:
raw[raw['heduc'].isnull()].shape
Out[53]:
In [64]:
raw = pd.read_csv("botswana.tsv", sep="\t", index_col=False)
raw['nevermarr'] = np.where(raw['agefm'].isnull(), 1, 0)
raw = raw.drop('evermarr', 1)
raw['agefm'].fillna(0, inplace=True)
raw.loc[(raw['nevermarr'] == 1) & (raw['heduc'].isnull()), 'heduc'] = -1
raw[raw['heduc'].isnull()].shape
Out[64]:
In [65]:
raw.head()
Out[65]:
In [66]:
raw['idlnchld_noans'] = np.where(raw['idlnchld'].isnull(), 1, 0)
raw['idlnchld'].fillna(-1, inplace=True)
raw['heduc_noans'] = np.where(raw['heduc'].isnull(), 1, 0)
raw['heduc'].fillna(-2, inplace=True)
raw['usemeth_noans'] = np.where(raw['usemeth'].isnull(), 1, 0)
raw['usemeth'].fillna(-1, inplace=True)
raw.shape
Out[66]:
In [69]:
raw.dropna(subset=['knowmeth', 'electric', 'radio', 'tv', 'bicycle'], how='any', inplace=True)
raw.shape[0]*raw.shape[1]
Out[69]:
In [71]:
Out[71]:
In [76]:
m1 = smf.ols('ceb ~ age + educ + religion + idlnchld + knowmeth + usemeth +'\
'agefm + heduc + urban + electric + radio + tv + bicycle + nevermarr + idlnchld_noans +'\
'heduc_noans + usemeth_noans',
data=raw)
fitted = m1.fit()
print fitted.summary()
In [78]:
fitted.model.exog
Out[78]:
In [80]:
print 'Breusch-Pagan test: p=%f' % sms.het_breushpagan(fitted.resid, fitted.model.exog)[1]
In [87]:
m2 = smf.ols('ceb ~ age + educ + religion + idlnchld + knowmeth + usemeth +'\
'agefm + heduc + urban + electric + radio + tv + bicycle + nevermarr + idlnchld_noans +'\
'heduc_noans + usemeth_noans',
data=raw)
fitted = m2.fit(cov_type='HC1')
print fitted.summary()
In [88]:
print 'Breusch-Pagan test: p=%f' % sms.het_breushpagan(fitted.resid, fitted.model.exog)[1]
In [89]:
m3 = smf.ols('ceb ~ age + educ + idlnchld + knowmeth + usemeth +'\
'agefm + heduc + urban + electric + bicycle + nevermarr + idlnchld_noans +'\
'heduc_noans + usemeth_noans',
data=raw)
fitted = m3.fit()
print fitted.summary()
In [90]:
print 'Breusch-Pagan test: p=%f' % sms.het_breushpagan(fitted.resid, fitted.model.exog)[1]
In [91]:
m3 = smf.ols('ceb ~ age + educ + idlnchld + knowmeth + usemeth +'\
'agefm + heduc + urban + electric + bicycle + nevermarr + idlnchld_noans +'\
'heduc_noans + usemeth_noans',
data=raw)
fitted = m3.fit(cov_type='HC1')
In [92]:
print "F=%f, p=%f, k1=%f" % m2.fit().compare_f_test(m3.fit())
In [94]:
m4 = smf.ols('ceb ~ age + educ + idlnchld + knowmeth + '\
'agefm + heduc + urban + electric + bicycle + nevermarr + idlnchld_noans +'\
'heduc_noans',
data=raw)
fitted = m4.fit(cov_type='HC1')
m3.fit().compare_f_test(m4.fit())
Out[94]:
In [95]:
m3 = smf.ols('ceb ~ age + educ + idlnchld + knowmeth + usemeth +'\
'agefm + heduc + urban + electric + bicycle + nevermarr + idlnchld_noans +'\
'heduc_noans + usemeth_noans',
data=raw)
fitted = m3.fit(cov_type='HC1')
print fitted.summary()
In [ ]: